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Abstract 

In this paper we investigate some practical aspects concerning the 
use of the Restricted-Denominator (RD) rational Arnoldi method for the 
computation of the core functions of exponential integrators for parabolic 
problems. We derive some useful a-posteriori bounds together with some 
hints for a suitable implementation inside the integrators. Numerical ex- 
periments arising from the discretization of sectorial operators are pre- 
sented. 

1 Introduction 

For the solution of large stiff problems of the type 

u'it)^ f{yit))^Lu{t)+N{u{t)), (1) 

where L e K^^-'^ arises from the discretization of unbounded sectorial opera- 
tors and TV is a nonlinear function, in recent years much work has been done 
on the construction of exponential integrators that might represent a promising 
alternative to classical solvers (see e.g. [23] or [TS] for a comprehensive survey). 
As well known the computation of the matrix exponential or related functions of 
matrices is at the core of this kind of integrators. The main idea is to damp the 
stiffness of the problem (assumed to be contained in L) on these computations 
so that the integrator can be explicit. 

Under the hypothesis that the functions of matrices involved are exactly 
evaluated, the linear stability can be trivially achieved for both Runge-Kutta 
and multistep based exponential integrators and hence highly accurate and sta- 
ble integrators can be constructed. On the other hand, the main problem with 
this class of integrators is just the efficient computation of such functions of 
matrices, so that, very few reliable codes have been written (we remember the 
Rosenbrock type exponential integrators presented in [3], [T7], [SO])- For this 
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reason many authors are still doubtful about the potential of exponential inte- 
grators with respect to classical implicit solvers even for semilinear problem of 
type ©. 

An exponential integrator requires at each time step the evaluation of a 
certain number (depending on the accuracy) of functions of matrices of the 
type ip^{hL)v, where 

^^{hX) = exp(/iA), (2) 
Vk+i{h\) = — for fc = 0,1,2,... , 

being h the time step. Actually this represents the general situation for the 
Exponential Time Differencing methods, that is, the methods based on the 
variation-of-constants formula; for Lawson's type method (also called Integrat- 
ing Factor methods) only the matrix exponential is involved. We refer again to 
[24] and the reference therein for a background. 

Among the existing techniques for the computation of functions of matrices 
(we quote here the recent book of Higham [15 for a survey), in this context 
the Restricted-Denominator (RD) Rational Arnoldi algorithm introduced inde- 
pendently in [37, and 261 for the computation of the matrix exponential seems 
to be an reliable approach. It is based on the use of the so called RD rational 
forms, studied in [29] for the exponential function, 

where qi is a polynomial of degree < i. We refer again to [26] for the basic 
references about the properties and the use of such rational forms. While in 
the matrix case, the use of these approximants requires the solution of lin- 
ear systems with the matrix (/ — (5L), as shown in [30] in the context of the 
solution of ll]) when L is sectorial so typically sparse and well structured this 
linear algebra drawback can be almost completely overtaken organizing suitably 
the step-size control strategy and exploiting the properties of the RD Arnoldi 
method concerning the choice of the parameter 5. In other words the number 
of linear systems to be solved can be drastically reduced with respect to the 
total number of computations of functions of matrices required by the integra- 
tor. Therefore the mesh independence property of the method, that leads to a 
very fast convergence with respect to a standard polynomial approach (see again 
[26] )■ can be fully exploited for the construction of competitive integrators. 

A problem still open is that inside the integrator the rational Arnoldi al- 
gorithm (responsible for most of the computational cost) have to be supported 
by a robust and sharp error estimator. In the self-adjoint case the problem 
has been treated in [27] where the author presents effective a-posteriori error 
estimates, even in absence of information on the location of the spectrum of L. 
Anyway, in the general case, when ([1]) arises for instance from the discretization 
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of parabolic problems with advection terms and/or non-zero boundary condi- 
tions the numerical range of L, that we denote by F{L), may not reduce to a 
line segment. In this sense the basic aim of this paper is to fill this gap providing 
error estimates for the non-symmetric case using as few as possible information 
about the location of F{L). It is necessary to keep in mind that a competitive 
code for ([!]) should also be able to update L (interpreted as the Jacobian of /, 
[35] ■ [4]) so that F{L) is may be not fixed during the integration, and so it is im- 
portant to reduce as much as possible any pre-processing technique to estimate 
F{L). In particular assuming that F{L) C C~we shall provide a-posteriori error 
estimates for the RD Arnoldi process using only information about the angle of 
the sector containing F{L), angle that is typically independent of the sharpness 
of the discretization and hence computable working in small dimension. 

The paper is organized as follows. In Section 2 we present the basic idea of 
the RD rational Arnoldi method and in Section 3 we derive some first general 
error bounds based on the standard approaches. In Section 4, exploiting the re- 
lation between the derivatives of the function e^^^ and the Laguerre polynomials 
extended to the complex plane, we derive some a-posteriori error bounds. The 
problem of defining reliable a-priori bounds is investigated in Section 5. Sec- 
tion 6 is devoted to the analysis of the generalized residual as error estimator, 
that can be used to obtain information about the choice of the parameter 6 for 
the rational approximation. In Section 7 we present some numerical examples 
arising from the discretization of a one-dimensional advection-difFusion model. 
In Section 8 we provide some hints about the use of the RD rational Arnoldi 
method inside an exponential integrator with the aim of reducing as much as 
possible the number of implicit computations of (/ — dL)~^. Finally, in Section 
9 we furnish a deeper analysis concerning the fast rate of convergence of the 
method, that will provide further information about the optimal choice of the 
parameter S. 

2 The RD rational Arnoldi method 

In what follows we denote by ||-|| the Euclidean vector norm and its induced 
matrix norm. As already mentioned, the notation F{L) indicates the numerical 
range of L, that is. 



while the spectrum of L is denoted by (t{L). The notation 11^ indicates the 
space of the algebraic polynomials of degree < m. 



be the unbounded sector of the left half complex plane, symmetric with respect 
to the real axis with vertex in and semiangle 6. Let moreover Tg be the 




Given < 6* < f , let 



Se = {X:\aTg{-X)\<e} CC 



(3) 
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boundary of Sg. Throughout the paper we assume that F{L) C int{Sg), the 
mterior of Sg- Accordmgly, L is a so-called sectorial operator (see e.g. [19] 
Chap. V, for a background). 

Given a vector v G M^^, with ||w|| = 1, consider the problem of computing 

y« = ^,{hL)v, (4) 

where iff. is defined by The RD rational approach seeks for approximations 
to ipf.{h\) of the type 

■Rm-i,m-i(A) = (X -"jA^)'"^^ ' Pfc,m_i(A) e n™_i, m > 1, 

where ^ > is a suitable parameter. Turning to the matrix case, y'^^^ is approx- 
imated by elements of the Krylov subspaces 

Km{Z, v) — span |w, Zv, Z^v, Z™^^v^ , m > 1, 
with respect to v and the matrix Z defined by the transform 

Z = {I-SL)-^. 

In this sense the idea is to use a polynomial method to compute y^'^^ = fk{Z)v, 
where 

is singular at 0. 

For the construction of the subspaces Km{Z,v) we employ the classical 
Arnoldi method. As is well known it generates an orthonormal sequence {vj }j>]^ , 
with vi = V, such that KmiZ, v) — span {vi, V2, Vm}- Moreover, for every m, 

(5) 

where Vm — [vi,V2, ...,u,„], H,n is upper Hessenberg matrix with entries hij = 
Zvj and Cj is the j-th vector of the canonical basis of K™. 
The m-th RD-rational Arnoldi approximation to y'^^^ is defined as (see [2D] ) 

yi^^ =VMHm)ei. (6) 

It can be seen that 

yln^ =Pk,rn-l{Z)v, (7) 

where pj. G Ilm^i interpolates, in the Hermite sense, the function fk{z) in 
the eigenvalues of H,n (see [32] ) . 

As mentioned in the Introduction this technique has been introduced inde- 
pendently in [57] and [25]. Anyway, the idea of using rational Krylov approx- 
imations to matrix functions was originally introduced in [5]. More recently 
this approach has been extended to the case of multiple poles and is commonly 
referred to as RKS (Rational Krylov Subspace) approximation (see [21], [31], 

my 
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3 General error bounds 

Before stating a general error bound for the method, we need to locate F{Z). 
Consider the function xW = (1 — <5A)~^. Denoting by -Di/2,1/2 the disk centered 
in 1/2 with radius 1/2, let 

Gg = {z:z = x{\),XeSe}CD,^2A/2- (8) 

Its boundary. Eg, is made by two circular arcs meeting with angle 29 at and 
1. Regarding the field of values of Z, F{Z), we can state the following result 
that will be used frequently throughout the paper. 

Proposition 1 If F{L) C intiSe) then F{Z) C int{Ge). 

Proof. Obviously a{Z) — x(o'(i)), so F{Z) cannot lie entirely outside 
Gg. Now assume that there exists A e such that x(A) € F{Z), that is, 
F{Z) n Eg 7^ 0. Hence, there exists y e C*^, \\y\\ = 1, such that 

yH{I-6L)-'y^^^. (9) 
Defining x := {I — SL)~^ y we easily obtain 

1 



x" (l-dL^)x^ 



and hence 



By ^ we have 



1-5X' 
^x^L^x 1 



5- 



(l-<5A)||x|r 



M\l-5X\>1. (10) 
Now let us define ^^^^^ e F{L). We have 



M'^{l-6\)-\l-5t,r\ (11) 



and hence 



that implies 



Im 



Im 



((1-<5A)-^(1-^m)"') -0, 
((1-5A)-^) _ Im((l-^M)"') 



Re((l-5A) Rc(^(l-5/i) 



-1 



(12) 



Now since (1 — (5/i) ^ G int{Gg) and (1 — 5X) ^ e Sg, by ([T^ it must be 
Re ((1 - (5A)"^) > Re ((1 - (5^*)"^) and Im ((1 - 5A)~^) > Im ((1 - (5/i)"^ 



so that 



\l-5ii\-^ < |l-,5ArV 
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Using this relation, by ((TT]) we finally have 

||a;||'|l-(5A|' < f, 

that contradicts ((TU]). Since the field of values is connected the proof is complete. 



Remark 2 In order to provide information about the geometry of F{Z), it is 
worth referring to \39f Theorem 5.2 in which the author proves that if L is an 
invertible matrix then 



lim 



1 



FiL). 



>oo VF((L-sJ)-i) 
Taking 6 = 1/s, we have that for small values of 6 

Going back to our method, the corresponding error Ek^m '■= y^^^ — Vm can 
be expressed and bounded in many ways (we quote here the recent papers [3] 
and [7] for a background on the error estimates for both polynomial and rational 
Arnoldi approximation to matrix functions). The following proposition states a 
general result. 

Proposition 3 Let G C 151/2,1/2 « compact such that F{Z) C int{G) and 
whose boundary T, is a rectifiable Jordan curve. For every p^^i £ Hm-i 



\Ek' 



< 



Ifkjz) -pm-i{z)\ 
2tt 7s dist{z,F{Z)) 



V - {zl - Z)VrnizI - H„i) ei \dz\. 

(13) 



Proof. Using the properties of the Arnoldi algorithm we know that for every 

Pm-l G Ilm-l, 

VmPm-liH,n)ei ^ Pm-l{Z)v . 

Hence from this identity it follows that, for to > 1 



Ek,m = fk{Z)v - Pm^l{Z)v - Vm{fk{Hra) - Pm-l{Hm))ei 



(14) 



Now since F{Hm) Q F{Z) we can write in the Dunford- Taylor integral 
form 



Eh 



1 



{fk{z) - Pm-l{z)) {zl - Z) ^ V -Vm{zl - Hrn) ^ El 



2TTi J 5] 

Collecting {zl — Z)^^ and using (see [33]) 

{zi-zy^ < 



dz. 



1 



dist{z,F{Z)y 
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we prove (fT^ . 
Now since 



{zI-Z)Vm(zI-Hmr^e^ 



qr,i{z)v 



where 

qm{z) = det(z/ - Hm), 

(see [25]), any bound for / |(7m(2)| and any choice for G and Pm-i 

leads to a bound for ||i?fc_m||- This technique has been used for instance in [25] 
and [16] . In particular in [26j the authors use the relation 

m 

hmiZ)v\\ ^Ylhj+ij, (15) 

J = l 

and the inequality 

\qrn{z)\>dtst{z,F{Z)r. (16) 

Going back to our situation, the main problem is that if we simply as- 
sume that F{L) C Sg (in other words F{L) arbitrarily large) we have that 
dist{z,F{Z)) — > as z — >■ (Re A — oo) because we have to consider the sin- 
gularity of fk at 0. Therefore using a lower bound like ([TC)) (but the situation 
remains true even for other approaches (cf. [IB])) terms of the type fk{z)/ z™'^^ 
would appear in (jl3p . In the exponential case (fc — 0) this is not a problem be- 
cause /o(2)/2:™"'"^ — > for z — ^ 0, but for fc > the situation changes completely 
since 

fk{z) S 1 

for z ^ 0. 

Because of the difficulties just explained, our approach for deriving error 
bounds is not based on the use of the Cauchy integral formula. Exploiting the 
interpolatory nature of the standard Arnoldi method, we notice, as pointed out 
also in [llj . that the error can be expressed in the form 

Ek^ni ^ gk,7n{Z)q,niZ)v, (17) 

where (cf. (0) 

, . fk{z) -Pk,m-liz) 

''^-^'^ ■= det(z/-g.„) ■ 

In [11] this relationship is used as the basis for the construction of restarted 
methods for the computation of matrix functions. 

We can state the following basic result that will be used throughout the 
paper and that allows to overcome the difficulties of working with formula (|T3)) . 

Proposition 4 Let F{L) C Sg and let t h/5. Then 



\Ek.m\\ < K 



T''{m + fc)! zeGe 



dz'' 



k{z)z 



(19) 



where 2 < K < 11.08. In the symmetric case we can take K ~ 1. 
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Proof. By [S] we know that 

\\9k.rn{Z)\\ < K max \gk,m{z)\ , 

zeF{Z) 

and hence by ([T5|) and pT|) 

nm 

Now, by induction one proves that for fc > 1 



hi^) - ^ ' (20) 



where 50(2) = 1 and 



T 



(z- 1)'= 



Sk{z) = Sfe_i(2;)z H e Ilfc for fc > 1. 

fc! 

Putting (gni) in (HI) we obtain 



gk,ra{z) 



t''{z- lY dei{zl - H^) 



Now, the polynomial t'^(z—1)''pj, ,„_j(z) e Ilm+fc-i interpolates in the Hermite 
sense the function fo{z)z'' — Sk-i{z)z in the eigenvalues of Hm and in z = 
1. Henceforth gk.m{z) is a divided difference that can be bounded using the 
Hermite-Genocchi formula (see e.g. [B]), so that 



|5fe,m(2)| < 



-'^(m + fc)! ^eco{{z,a{H^),l}) 



rn-\-k ' 



where co({z, cr (_//,„), 1} denotes the convex hull of the point set given by z, 
(j{Hm) and 1. Since a{H,n) C F{Z), and F{Z) C Ge by Proposition [U the 
result follows. ■ 



4 A posteriori error estimates 

By (|19p . in order to provide a-posteriori error estimates we just need to study 
the derivatives of the function /o(z)z'''. We need to introduce the generalized 
Laguerre polynomials, defined by 

4"'w-E(-i)<" + ")4- 

We can state the following result. 
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Lemma 5 Let t = ^ • For to > 1 
1 



d ■ 



Z ' 



r'=(TO + fc)! dz'"+*='^^ ;2m+fc+i -"'y ^(rn + k)l " ^^-^^ ^ ^ 



Proof. First of all remember that fo{z) = e^e Defining w = z/t and 
using Rodrigues' formula for Laguerre polynomials (see [T] p. 101) we obtain 



expf = 



2^ ^m+fc _^ 

exp (— cj^"^) (w^"^) 



1 



The result arises from the relation (see p. 240) 



L 



i-i-k) 

m-\-k 



(_) = (_i)-+i(_) 



fc+V'T^fc+l ("^- 1)! 



z (m + fc)! 



7- y-i-'-j ( ■ \ 



Before stating the main result we need to remember the following properties 
of the generalized Laguerre polynomials, that can be found in pp. 785-786. 



LI 



L2 



L3 



L(")(ziz,) = f](" + ")L5")(zi)^^(l-^2)' 



J=0 



exp( — ) L\'^\x) 



r(n + Q; + l) 
< — for X > 0. 



7i!r(a + l) 

Proposition 6 Given r > 0, let z = {l + Sre'^^) ^ G Eg. Let moreover 
Cj(6») (l + V2(l-cos6l)y . 



Then 



m — l 



^L^^I^t) < 4t^,_,.(r) c,(^), 



3=0 



T + hr v-^m-l /to + fc - 7 — 1\ 



(22) 

(23) 

(24) 
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Proof. For z = (l + Sre'") ^ 



-^T + hre'\ r>0. 
z 

Using LI with a = k, (3 — 0, zi — t and Z2 — hre^^ , and then L2 with zi — hr 
and Z2 = e**, we have 



m — 1 

J=0 



< E|^™Vi(-)|Eki°'(H 

j=0 



(25) 



Since 



E 



s=0 



^0,(9), 



formulas and ([M)) are obtained applying L3 to lI^' (/ir) and then to -^m-i-i (''') 



Theorem 7 Assume that F{L) C Sg, wit/i ^' < f • T'/ien 



II ^fe,^ 



< K- 



„T cos 6—m — k—l 



< K- 



j-ni-{-k 



2 cose - 1 

m+k+l 



2(m + A: + 1) 
2 cos 6* - 1 



1=1 

C^,™Wn'^'+i.- (27) 



where 



(m- 1)! 
(to + fc) 



J=0 



(to — 1)! ■r-^m-i /m + fc — J — 1 



1 



(to + fc)! '^j= 
and K defined as in Proposition 
Proof. For z e Se 



1 + 5re^\ r > 0, 



and 



/o(2;) = e"' ^ = e 



Hence, using (|T9|, ([21]) and ([23]) we obtain 

\\Ek^\\ < Xmax|e-'"^(™^''-5)(l + ^re''')"+'+'| 

r>0 I \ / I 



xi^;^|4^L,_,(r)|c,(e)n/^m.. 

i=l 



(to + fc)! 

^ ^ i=o 



(28) 

(29) 



(30) 
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Since for 9 < tt/3 



(looking for the maximum with respect to r), we immediately obtain (1261) . Using 
again (|T9l) and ([2T|) but now with p4|) we arrive at the coarser bound ((27l) . ■ 



Remark 8 While formulas \2b]) and \2T^ theoretically hold for < ^ since 
hm+i,m = for m < M , it is necessary to point out that for 6 k ^ we may 
observe a rapid growth of the term 



2 cose - 1 



Uh 

1=1 



depending of course on the problem, so that the bounds may be useless. This 
situation is caused by the hound i23\} that leads to the appearance of the term 
2 cos 9 ^ 1 at the denominator. Working in inexact arithmetics the situation is 
even more difficult because of the loss of orthogonality of the vectors Vj of the 
Arnoldi algorithm and hence the accumulation of errors on the entries hi+i^i. 
For these reasons, in practice, formulas \20jl and {21^ should he used only for 9 
not much close to ^. 

Remark 9 For the exponential case (k = 0) we have 
and hence by ^2j 



m+1 



1=1 



Remark 10 In the self-adjoint case 9 = we have Cj{9) — 1 and formula \21\j 
simplifies to 

\\Ek,m\\ < {k + iy. 11^^+1.- 

The reason for which we consider two bounds in Theorem[7]is that the second 
one (P7)) allows us to define suitably the parameter t (and then S) while the 
first one (j26p should be used whenever t has been defined. Indeed, assuming 
YliLi hi+i.i independent of S and then of r (actually this is not true as we 
explain in Section 9) by ([27l) . looking for the minimum of e'^'=°''^T~('"+*^^ we 
easily find that the optimal value for r is given by 

r-^. (32) 
cos 9 
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The a-posteriori bounds provided by Theorem [7] depends substantiahy on 
the semiangle of the sector containing F{L). Therefore, the most natural 
way to proceed is to compute the boundary of F{L) using the standard codes 
available in hterature (as for instance the Matlab code fv.m by Higham [H]). 
It is important to observe that 9 is generally independent of the discretization 
so that one can work in smaller dimension. 



While the hypothesis F{L) C Sg of Theorem [7] is extremely general, the 
underlying assumption is that L represents an arbitrary sharp discretization 
of an unbounded operator. On the other side, if it is known that F{L) is 
contained in a bounded sector then Proposition [3] can be used to derive sharper 
error estimates. In general we may refer again to [3| and the references therein 
for a background on the most used techniques based on the use of the integral 
representation of the error. 

Anyway, here we want also to show how to adapt our approach in precence 
of more information on F{L). Let Dq h be the disk centered at with radius 
R, and assume that F{L) C SgCi Dqh. Using again and (PTI) . we arrive at 
the bound 



\Ek,m\\ < K max 

0<s<hR 



(^-1)! 
(m + fc)! 



(33) 



In order to define a suitable value for r, we just need to bound the Laguerre 
polynomials as in (P^ . so that the optimal value is obtained ooking for the 
minimum of 

r(l + -) 

A good approximation for this minimum is given by 



r = y/2hR{m + k + l), 
that is obtained considering the bound 



(34) 



m+fe+l / hR 

< exp (m + k + 1) 



Using this value of t we can derive practical error bounds seeking for the maxi- 



mum of the function 
[0,hR]. 



(1 + f )"+'+' L$°)(se'^) (cf. 



25t) in the interval 



5 A-priori error bounds 

Formula (I32p obviously requires to know the number of iterations that are nec- 
essary to achieve a certain accuracy. In this sense we need to bound in some 
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way YlTLi By (fT31) and since 

1 1 Qm 

{Z)v\\ < \\p,n{Z)v\\ 

for each monic polynomial Pm of exact degree m (see [35] P- 269), a bound for 
nl^li f^i+i.i can be stated using Faber polynomials as explained in [2 , that leads 
to 

nm 
= ||g„.(Z)t;|| < 27(0)", (35) 
2 — 1 

where 7(0) is the logarithmic capacity of a compact G containing F{Z) and 
where fk is analytic. 

Proposition 11 Let 9* = 0.48124 and assume that F{L) C Se, with 9 < 9* . 
Then for t = (to + fc)/ cos 9 

\\Ek,^\\ < nKp{9r, (36) 

where 

p(g):=(l + V2(l-cosg)) ^ <1 V 0<9<9*. (37) 

\ / 4 cos 9 — 2 TT — 9 

Proof. Since F{Z) C Gg by Proposition [TJ let us consider the compact 
subset G — Gg. The associated conformal mapping 



is given by 



V' : C\{w : \w\ < 1} ^ C\Gg, 



(w + 1)^ 



2-1/ ■ 



{w + 1)2-" - {w - 1) 

1 1 1 (1 - i^) (3- i^) 1 ^ / 1 , 



2(2 -J/) 2 6 2 ~iy w 



where v = 29 /tt. The coefficient of the leading term of the Laurent expansion 
(|38|) is the logarithmic capacity, so that by (1351) we have 



/ 1 \ ™ 

n;:/-M<2U7^ . (39) 



.2(2 - ly) 

Inserting this bound in (P7)) we easily obtain for 6* < |- 

m + A: + 1 / cos ^ 



where the second inequality arises from the choice r = (to + k)/ cos 9. 
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Now, by the definition ([2^ . it is rather easy to show that 
(m — 1)! fm + k — j — 1 



1 ( m-1 Y 
- k\{m + k)\m + k-l) 



so that 



/c+l 

cost/ \ k+3 



' " - A:!cos6l V 2cos6'- 1 ' ^ ' 



Since 1/2 <$ (61) < 1 for < 6* < 6**, and since for each fc > 

( COS0 y''\k+3 < ^ cosr 
A:!cos6' l^2cos6'- ly " cos6i* V2 cos6»* - 1 

the proof is complete. ■ 

Remark 12 Proposition \11\ shows the mesh-independence of the method for 
9 < 9* since the bound i36\) is independent of the discretization of the underlying 
sectorial operator. In Section 9 this considaration is extended to 9 < ^ . By ^Jd^ 
and Ji^Tp , in the self-adjoint case (K — 1) the bound h36\) reads 

It is worth noting that by (|14p for every p„i-i G ^m-i we have that 
||-Bfe,„J < 2i^max|/fc(z) -p„i-i(z)| , 

where we assume that G C -Di/2,1/2 is compact, connected, with associated con- 
formal mapping 0, and such that F{Z) C G. Therefore, in principle, one could 
try to derive a-priori error bounds choosing suitably the polynomial sequence 
{Pm-i\ra>i- ^^lyway, the classical results in complex polynomial approximation 
state that even taking {Pm-i\ rn>i ^ sequence of polynomials that asymptot- 
ically behaves as the sequence of polynomial of best uniform approximation of 
fk on G (see e.g [33] for a theoretical background and examples) we have 



max|/fc(z) -pm-i{z)\ 



1/m ^ 



— as m cx), 
R 



where R is such that </>(— i?) 0, since fk is singular at {maximal convergence 
property, see e.g [38] Chapter IV). The main problem is that assuming L to be 
unbounded, G G and consequently R = 1. 

For this reasons, in our opinion the only reasonable approach to derive a- 
priori error bounds, is to define {pm-i}m>i as a sequence of polynomials in- 
terpolating fk at point belonging to G, and then to use the Hermitc-Genocchi 
formula to bound the divided differences. Using this formula and taking for 
instance Pm-i as the sequence of interpolants at the zeros of Faber polynomials 
we just obtain the error bound given in Proposition [TT] (see [25)). 
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6 The generalized residual 

By the integral representation of function of matrices and ([6]), we know that the 
error can be written as 

Ek,m = :^f fk{z)[{zI-Z)-'v-Vm{zI-H^)-'ei]dz. (41) 

In order to monitor the approximations during the computation we can consider 
the so-called generalized residual [I^, defined as 

Rk,m = ^ J fk{z)r.m{z)dz, (42) 

which is obtained from (j4ip by replacing the error 

[zl - Z)-^v -V,r.{zl - H^)-^e^ 
with the corresponding residual 

r„(z) = v - {zl - Z)Vrn{zI - F„)"^ei. 
Using the fundamental relation ([5]) we have immediately 

r„(z) = h„i+i,m{e"{zl - Hrn)~^ei)vm+ij 
and inserting this relation in (|42|) we obtain 

SO that we may assume 

Ek,m ~ |l-Rfe,m|l = h,n+l,m | e^/fe (i?m)ei | . (43) 

In order to show the reliability of this approximation let us consider the 
operator 

Lu = -u" + cu', c>0, (44) 

discretized with central differences in [0,1] with uniform mesh h — 1/(M + 
1), and Dirichelet boundary conditions. For our examples, we consider the 
computation of ipj.{hL)v for fc = 1,2, with v = (1, Vf^ /^/M, comparing the 
exact error and the generalized residual. We take Af = 1000, h — 0.1, and we 
consider the cases of c = 2 and c = 4, whose corresponding sector semiangles 
are respectively 9 = 0.201 and 6 — 0.425. We define r = 15/ cosO. The results, 
collected in Figure 1, shows the accuracy of the approximation 

It is necessary to point out that the use of (|43p has the basic disadvantage 
that it requires the computation of fk{Hm), "ni ~ 1,2,..., and this is a com- 
putational drawback whenever a great amount of matrix functions evaluations 
are required to integrate a certain problem, even if m can be considered much 
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smaller than M . Moreover, it frequently happens (as in our experiments) that 
the generalized residual tends to underestimate the error during the first itera- 
tions, and this can be particularly dangerous when computing Lpf.j^-^{hL)v with 
1 1 11 1 1 <C 1, as for instance in the case of the computation of the internal stages of 
an exponential Runge-Kutta method, in which \\v\\ = 0{h). 

On the other side, exploiting the mesh independence of the method the 
generalized residual can be successfully used to estimate the optimal value for 
the parameter t, that is Topt = {m-\- k) / cos9. In other words, using a coarser 
discretization of the operator we look for the value of m such that using the 
corresponding Topt we obtain a certain tolerance in exactly m iterations. For 
the experiments reported in Figure 1 we considered the discretization of (|44p 
with only M = 50 internal points, observing in both cases that |i?m| < le — 12 
for 771 > 13. For this reason we have chosen t — 15/ cos9. 



Figure 1 - Comparison between the exact error and the generalized residual for problem (|44|) with h — 0.1. In 

7 Numerical experiments for the a-posteriori er- 
ror bound 

For our numerical experiments we consider again the operator (I44p . discretized 
as in previous section. We consider the computation of the functions ipi.{hL)v, 
with V as before and k = 0,1,2, for h = 0.5 (Figure 2) and h = 0.05 (Figure 
3). In all examples we do not consider the symmetric case corresponding to 
c = (already investigated in ^7\), but only the cases c = 2 and c — A. As 
before, for the choice of r we examined the behavior of the method for the 
coarser discretization of the same operator with only M = 50 interior points, 
thus exploiting the mesh independence of the method. The analysis suggested 
to take T — 8/ cos 6* for all experiments with h ~ 0.5 and r — 15/ cos 6* for those 
with h = 0.05, thus independently of the function and c, using the tolerance 
le — 12. Inside the Arnoldi iterations the vectors Zvj, j > 1 (cf. Section 2), are 
computed via the LU factorization of / — 5L. The error bound is given by (1^ . 
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Figure 2 - Error and error bound p6| for = 0, 1, 2, h = 0.5, L arising from (|44l) with c = 2 and 



17 



Figure 3 - Error and error bound ([26l) for A: = 0, 1, 2, h — 0.05, L arising from (j44|) with c = 2 and c = 4. 

Comparing Figure 2 with Figure 3 we can observe that the method tends to 
become sfower reducing h. The reason is that for small values of h, the rate of 
the decay of the singular values of Z becomes slower and this reduces the rate 
of the decay of Hi" i hi+n. A deeper analysis of this behavior will be presented 
in Section 9. 

8 Non-optimal choice of r 

Employing the RD Arnoldi method inside an exponential integrator requires 
some considerations. First of all, in our opinion the method can be used only 
if the implicit computation of Z can be obtained with a sparse factorization 
technique. The use of an inner-outer iteration can be too much expensive in 
this context. Indeed, the basic point is that organizing suitably the code one 
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can heavily reduce the number of factorizations of I — SL (see e.g [30]), because 
the method seems to be really robust with respect to the choice of t. For this 
reason we want here to show what happens taking r even quite far from the 
optimal one. 

For simplicity (the situation is representative of what happens in general) let 
us assume to work with exponential function and = 0. We assume moreover 
that the corresponding bound (|3ip (in which Cj{9) = 1, j > 0) is equal to a 
prescribed tolerance for a certain m with the theoretical optimal choice Topt = 
m. We seek for the interval Im,n = [ti,T2] such that for t e the number 
of iterations necessary to achieve the same tolerance is at most n {> m). Using 
(|5T|) and the approximation hm+i.m « 1/4 (to > 1) that is obtained forcing the 
equal sign in the a-priori bound (|39p . in Figure 4 we can observe the result for 
n = TO + 1, m + 2. For each to the corresponding extremal points ti and T2 of 
the intervals Im,m+i and Im.m+2 are plotted. These points are obtained solving 
with respect to r the equation (cf. (PTjl ) 

^^(2(n+i)r+in. /.+M = — (2(m + i)r+in. 

T -'-J-i=l TO -1-^-1=1 

for n = TO + 1, TO + 2. 



Figure 4 - Boundary of the region Im,m+i and Im,m+2- 

We point out that the results are even a bit conservative with respect to 
what happens in practice, and this is due to the approximation hm+i,m ~ 1/4. 
Indeed larger intervals would be obtained taking hm+i,m < 1/4 as it occurs in 
practice. 
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In order prove the effectiveness of the above considerations let us consider 
again the operator (l44l) with the usual discretization. We consider the case 
c — 2, k — I for h — 0.1. To define r we consider again the discretization 
with M = 50 interior points observing the generalized residual. This leads us 
to define r = (m + k)/ cosO with m = 14. In Figure 5 we consider the behavior 
of the method for r, r/2 and 2r. 



Figure 5 - Error and error bound (pS)) for fc = 1, /i = 0.1 and L arising from (^H) with c — 2. Method applied v 
T = 15/ cose*, r/2 and 2t. 

The robustness of the method with respect to the choice of r is maybe the 
most important aspect concerning its use inside an exponential integrator. We 
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want to give here some practical suggestions assuming to use a sparse factoriza- 
tion technique to solve the Unear systems with / — 6L, that, computationally, 
has to be considered the heaviest part of the method. 

1. Working in much smaller dimension compute 9 and use the generalized 
residual to estimate the initial Topt- 

2. For nonlinear problems, interpreting L as the Jacobian of the system ([4], 
[35j). it is necessary to introduce some strategies in order to reduce as much 
as possible the number of updates of L during the integration, since each 
update would also require to update the factorization. As for exponential 
W- method (see [17], [30]), we suggest, whenever it is possible, to work 
with a time-lagged Jacobian and hence to introduce the necessary order 
conditions in order to preserve the theoretical order. 

3. Using a quasi-constant step-size strategy (without Jacobian update) allows 
to keep the factorization oi I — 5L constant for a certain number of steps. 
Whenever it is necessary to update the stepsize hoid — > hnew without 
changing the Jacobian, if we want to keep the previous factorization of 
I — SoidL we just need to consider the ratio r = hnew / Sold- If (indicatively) 
it is bigger than 2Topt or smaller than Topt/ '2' (cf. Figure 4 and 5), where 
Topt arises from a previous analysis of the generalized residual, then we 
need to update the factorization (cf. again [30]), otherwise we can keep 
the previous one. In this phase, however, one can even considers other 
strategies to define suitably the window of admissible values of t around 
Topt, taking into account of the local accuracy required by the integrator, 
the norm of v, etc. 



Looking carefully at Figure 5 we notice that while the analysis in smaller dimen- 
sion suggested to take t = 15/ cos 9 for reaching the desired tolerance in exactly 
14 iterations the method is unexpectedly a bit faster taking n = t/2 (sec- 
ond picture). The analysis was correct because in larger dimension the method 
actually achieves the tolerance in 14 iterations (first picture). In order to under- 
stand the reason of this behavior, we need to remember that the definition of 
Topt = {m + k)/ cos 9 given at the end of Section 4 was based on the assumption 
that YiiLi is independent of i5, but this is not true. In what follows we 

try to provide a more accurate analysis studying the decay of Hl^i 

We denote by (Tj, j > 1, the singular values of Z. Moreover we denote by 
Aj, j > 1 the eigenvalues of Z and assume that |Aj| > |Aj+i| for j > 1. We have 
the following result (cf. [28] Theorem 5.8.10). 

Theorem 13 Assume that 1 ^ cr(Z) and 
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The superlinear decay of H^li ^i+i,i 




(45) 
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iet = nlli(^ ~ ^0- Then 

i/p 

m 

where 



,AZ)\\<{^T", (46) 



As already shown in Section 4 

nm 
. , hi+i,^ < \\p„i{Z)v\\ 

for each monic polynomial pm of exact degree m (see p. 269), so that The- 
orem [13] reveals that the rate of decay of YYiLi ^i+i.j is superlinear and depends 
on the p-summability of the singular values of Z. We remark moreover that an 
almost equal bound has been obtained in [13] studying the convergence of the 
smallest Ritz value of the Lanczos process for self-adjoint compact operators. 

In practice, the use of (|46|) requires the knowledge oip and a bound for 77, that 
is, information about the singular values of the operator Z. As a model problem 
we consider again the operator L defined by (H^ with c — 0, whose eigenvalues 
are (jtt)^, j > 1, so that the eigenvalues of Z are given by Xj — 1/(1 -I- S (jtt)^). 
In this case (|45l) holds for l/2<p<lso that Z can be referred to as a trace 
class operator (see again [28]). Hence, taking for instance p — 1 we have 



E-r < 

i>i 




arctan ( v^tt 



1 f 1 

■ arctan f 



and so 



n: 



< (47) 

h^+l,^ < \\Prn{Z)v\\ < ( ^) . (48) 



m 



The bound (|48|) reveals that the rate of decay depends on the choice of 5 and 
then on h. For large values of 5, say 5 > 1, the bound (|47|) can be heavily 
improved exploiting the properties of the arctan function and the convergence 
is extremely fast. The following proposition states a general superlinear bound 
that can be used when L is an elliptic differential operator of the second order, 
so with singular values growing like p . The proof is straightforward since we 
just require to bound X]j>i '^^ji ^^^d apply P5|) with p = I. 

Proposition 14 Let L be an elliptic differential operator of the second order. 
Then there exists a constant C such that 

{ C \ 

< . (49) 



n: 



m 
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This proposition can easily be generalized to operator of order s > 1, ex- 
ploiting [28j Corollary 5.8.12 in which the author extends Theorem [T3l for p > 1. 
Anyway, this is beyond the purpose of this section. 

From a practical point of view, formula (|49p is almost useless since too much 
information on L would be required. On the other side, it is fundamental to 
understand the dependence on 5. Setting as usual t — h/5 and putting the 
corresponding bound (|49)) in Theorem [7] (formula ((27|) '). we easily find that the 
theoretical optimal value for r is obtained seeking for the minimum of 

gT COS 6 — — k — 1 

with respect to t, that is, 

m + 2fc 

This new value, less than (m + fc) / cos0, explains our considerations about 
Figure 5 given at the beginning of this section. 

We need to point out that since the choice of Topt is independent of C and 
/i, formula (|49|) is quite coarse for small values of h and not able to catch the 
fast decay of JlI^Li I^i ^-^y case, if an estimate of C is available an a-priori 

bound for the RD Arnoldi method can be obtained taking 

n:.'-— ■"{(7£)'"'<2(^)"'}' 

(cf. (|39I) V Consequently we argue that 

m + 2k m + k 

^ ^ opt 1 

2 cos 9 cos 

with Topt close to (m + 2k) / (2cos0) for h large and to (m + k) / cos 6 for h 
small. 

10 Conclusions 

In this paper we have tried to provide all the necessary information to em- 
ploy the RD Arnoldi method as a tool for solving parabolic problems with 
exponential integrators. The little number of codes available in literature, and 
consequently, the little number of comparisons with classical solvers is a source 
of skepticism about the practical usefulness of this kind of integrators. Indeed, 
with respect to the most powerful classical methods for stiff problems, the com- 
putation of a large number of matrix functions, generally performed with a 
polynomial method, is still representing a drawback because of the computa- 
tional cost. The use of polynomial methods for these computations may even be 
considered inadequate whenever we assume to work with an arbitrarily sharp 
discretization of the operator, since this would result in a problem of polyno- 
mial approximation in arbitrarily large domains. For these reasons, the use of 
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rational approximations as the one here presented, should be considered a valid 
alternative because of the fast rate of convergence and the mesh independence 
property, provided that we are able to exploit suitably the robustness of the 
method with respect to the choice of the poles, as explained in Section 8 for our 
case 
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